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Abstract 



Recent experimental and theoretical results have shown the existence of a 
liquid- liquid phase transition in isotropic systems, such as biological solutions 
and colloids, whose interaction can be represented via an effective potential 
with a repulsive soft-core and an attractive part. We investigate how the 
phase diagram of a schematic general isotropic system, interacting via a soft- 
core squared attractive potential, changes by varying the parameters of the 
potential. It has been shown that this potential has a phase diagram with 
a liquid-liquid phase transition in addition to the standard gas-liquid phase 
transition and that, for a short-range soft-core, the phase diagram resulting 
from molecular dynamics simulations can be interpreted through a modified 
van der Waals equation. Here we consider the case of soft-core ranges compa- 
rable with or larger than the hard-core diameter. Because an analysis using 
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molecular dynamics simulations of such systems or potentials is too time- 
demanding, we adopt an integral equation approach in the hypernetted-chain 
approximation. Thus we can estimate how the temperature and density of 
both critical points depend on the potential's parameters for large soft-core 
ranges. The present results confirm and extend our previous analysis, showing 
that this potential has two fluid-fluid critical points that are well separated in 
temperature and in density only if there is a balance between the attractive 
and repulsive part of the potential. We find that for large soft-core ranges 
our results satisfy a simple relation between the potential's parameters. 

61.20.Gy, 61.20.Ne, 61.25.Mv, 64.60.Kw, 64.70.Ja, 64.60.My 



2 



Typeset using REVT^ 



I. INTRODUCTION 



The phase diagram of a typical monatomic substance is comprised of sohd and fluid 
phases, with the fluid phase separating below the critical point into gas and liquid phases. 
The prototype of such substances are simple (i.e., argon-like) fluids. Interparticle inter- 
actions in these systems can be appropriately described by the well-known Lennard- Jones 
potential. Other simple models — such as those described by the hard-sphere square-well po- 
tential or by the hard-sphere- Yukawa potential — exhibit similar phase diagrams. All these 
potentials consist in a short-range harshly repulsive core plus a longer-ranged attraction. 
New insights into the relationship between phase diagrams and interparticle interaction 
emerged recently from the finding that when the range of the attractive component is suf- 
ficiently small, the liquid phase and the gas-liquid critical point become metastable with 
respect to crystallization^"^. Shouldered potentials, i.e., potentials with a hard core and a 
finite repulsive shoulder, exhibit even more exotic phase diagrams. Simulations and theories 
showed that such potentials may give origin to non-trivial phase behaviors, such as isostruc- 
tural sohd-sohd transitions and liquid-hquid transitions^"^^. The key to this complex phase 
behavior resides in the peculiar penetrability of the repulsive core, a feature that gives rise 
to a density- dependent effective interaction. 

The possible existence of a liquid-liquid phase transition for single-component systems 
in particular has received considerable attention in recent years. A direct evidence of this 
phenomenon has been observed on the experimental side in liquid phosphorous^^'^^. Exper- 
imental data consistent with a liquid-liquid phase transition have been presented for other 
single-component systems such as water^^^^°, silica^^'^^, carbon^^, selenium^"^, and cobalt^^, 
among others. A recent theory has shown that a phase transition occurring in a solution 
of DNA-coated colloids can be understood in terms of the liquid-liquid phase transition^^. 
A liquid-liquid critical point has also been predicted by simulations for specific models of 
supercooled water^'''"^^, carbon^°, phosphorous^^, supercooled silica^^'^^'^^, and hydrogen^^. 

We have recently shown through molecular dynamics (MD) simulations^^ that a system 
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of particles interacting through an isotropic potential with an attractive well and a repulsive 
component consisting of a hard core plus a finite shoulder may possess a high-density liquid 
phase and a low-density liquid phase. Potentials with such characteristics were used to model 
interactions in a variety of systems including liquid metals, metallic mixtures, electrolytes, 
and colloids, as well as anomalous hquids, hke water and sihca^^"^^. 

In spite of the simplicity of the model, the physical mechanism that causes the liquid- 
liquid transition for a potential with a hard core plus a repulsive shoulder and an attractive 
well is not easy to assess since it arises from an interplay of the different components of 
the interparticle interaction. To disentangle the role of each component it is necessary to 
investigate the dependence of the phase stability of the fluid on the potential parameters. 
This task was undertaken in Ref.^^, where the results of MD calculations performed for 
several sets of parameters were presented. The resulting behavior of the critical points was 
interpreted through a modified van der Waals equation (MVDWE)^^, a mean field approach 
assuming that the effect of the repulsive shoulder at different densities p and temperatures 
T can be taken into account by an effective excluded volume depending on both p and T. 
However, the analysis was hmited to cases where both the soft-core range and the attractive 
range are smaller than the hard-core range a and the total interaction range does not exceed 
2.6a. Indeed, MD becomes quite time-demanding for larger interaction ranges. Nevertheless, 
there are cases such as biological solutions and colloids where the soft-core range could be as 
large as the hard-core or even larger^^. For this reason and to gain a better understanding of 
the role played by each component of a soft-core's attractive potentials, we will now explore 
how the phase diagram changes when the soft-core range exceeds the hard-core diameter. 

For this we use the hypernetted-chain (HNC) integral equation for the radial distribution 
function^^. This approach can be considered in many respects as an intermediate between 
MD and MVDWE. In fact, the HNC equation is far less time-consuming than a simulation, 
but is by no means as accurate. On the other hand, being a microscopic theory, it is not 
based on the assumption typical among mean field approaches such as the MVDWE that 
particles experience a uniform attractive potential. Hence it is intrinsically more accurate. 
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Since the HNC equation provides a fast way of estimating the position of the critical points, 
we perform an extensive investigation in the space of the potential parameters, considering 
an extremely ample number of combinations. Thus we are able to frame previous results into 
a wider perspective and obtain a better understanding of the physical mechanism leading 
to a liquid-liquid transition in one-component fluids. 

Our results show that the high-density critical point can be found only when there is 
a balance between the attractive part and the repulsive part of the potential. In Ref.^^ 
this balance was expressed through the mean fleld strength of attraction, a parameter pro- 
portional to the attractive range w^/o (Fig- 1) and inversely proportional to the repulsive 
energy Ur, for fixed attractive energy C/4. Here we find an approximated relation between 
Ur/Ua and wa/wr that is well- verified for large soft-core ranges and quantifies the ideal 
balance between the repulsive and the attractive components of the potential more efi^ec- 
tively. Our results show that the liquid-liquid phase transition could be found in systems 
with small repulsion if the attraction is small as well, with Ur/Ua wa/wr, and in systems 
with wide repulsion, with Ur/Ua oc 3wa/wr. Typical systems with these characteristics 
could be colloids, where the effective repulsion and attraction can be regulated^^'^^. 



II. THE ATTRACTIVE SOFT-CORE POTENTIAL 

A soft-core potential with an attractive interaction at large distances was first proposed 
by Hemmer and StelP^ to understand the possibility of the solid-solid critical point in ma- 
terials such as Ce and Cs and was studied through an exact analysis in ID. Other soft-core 
potentials with an attractive well were proposed and studied with approximate methods 
or with numerical simulations in 2D to rationalize the properties of liquid metals, alloys, 
electrolytes, colloids, and the water anomahes^^"^^'-'^^'^^'^^'^^'^^. 

The peculiarity of such potentials is the presence of two repulsive length scales. This 
feature is typical of systems with core-corona architecture such as, for example, star poly- 
mers. However, isotropic soft-core potentials have also been proposed as effective potentials 
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resulting from an average over the angular degrees of freedom for systems where the distance 
of the minimum approach between particles depends on their relative orientation. Thus, in 
some respects, they have been considered"^^'"^^'"^"^ as simplified models of complex anisotropic 
interactions, such as those resulting from the hydrogen bonding between water molecules. 

The model potential considered in this paper is similar to that investigated in 
Refs.^^'^^'^"^'^^ and is an isotropic pair potential with two characteristic short-range repul- 
sive distances: one associated with the hard-core exclusion between two particles and the 
second with a weak repulsion (soft-core), which can be overcome at large pressure. More 
precisely, our pair potential U (r) (Fig. 1) consists of a hard-core of radius a, a repulsive 
square shoulder of height Ur extending from r — a to r — b, and an attractive component 
having the form of a square well of energy —Ua<0 extending from r — bto r — c (here r is 
the interparticle distance). Choosing a and Ua as length and energy units respectively, this 
potential depends on three free parameters: the width of the soft-core WR/a = {b — a)/a, 
the width of the attractive well wa/ci = (c — b)/a, and the soft-core energy Ur/Ua- 

Our aim is to understand how the position of the critical points in the thermodynamic 
plane changes upon varying the parameter values. In Ref.^^ we investigated a number 
of cases with wr < a and wa < a through MD simulations and presented a mean field 
approach with an MVDWE to rationalize the results. However, the MD analysis of potentials 
with Wr > a would require very large computation times, therefore we study this case 
with integral equations in the HNC approximation, which represents a compromise between 
accuracy and economy. 

III. THE HYPERNETTED CHAIN INTEGRAL EQUATION APPROACH 

The spatial distribution of a system of particles may be conveniently described by the 
radial distribution function g{r)'^^, a quantity directly measurable by scattering experiments 
and related to the thermodynamic properties of the fluid. One of the theoretical approaches 
most used to calculate this function is represented by integral equations. These are based on 
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the so-called Ornstein-Zernike (OZ) relation which relates the total pair correlation function 
h{r) = g{r) — 1 to the direct correlation function c(r), which describes the contribution 
coming from the direct interaction between two particles at distance r. Both h{r) and c(r) 
are unknown functions, thus to solve the OZ relation, one needs another relation (closure) 
between these two functions. This is by necessity an approximate relation and it is obtained 
in the HNC equation^^ by neglecting the sum over a specific class of diagrams in the dia- 
grammatic expansion of g'(r)^^. In principle, this approximation is expected to work better 
at lower p, where the direct correlation function c(r) is more relevant than the correlation 
propagated through the other particles. 

The solution of the system formed by the OZ relation plus the HNC closure is obtained 
through a numerical iterative procedure which is stopped when the difference between two 
consecutive elements of the succession generated (in the space of the distribution functions) 
is smaller than some given value. Based on the standard iterative procedure, different al- 
gorithms can be used to improve the accuracy and rapidity of convergence of the HNC 
equation's numerical solution. However, independent of the algorithm used, there is a re- 
gion in the p-T plane where no solution can be found, i.e., for any p, there is a T below 
which the numerical algorithm does not converge. This defines an instability line (IL) of the 
theory in the p-T plane. The nature of the locus of instabilities of the HNC equation and 
its relationship with the spinodal line of the fluid was the object of extensive investigation^^ 
which showed that though the isothermal compressibility increases considerably when the 
instability line is approached from above, there is no real divergence of this quantity. Thus, 
identifying the IL of the HNC equation with the spinodal line of the fluid, which is charac- 
terized by a diverging compressibility, is not possible. Another well-known deficiency of the 
HNC is its thermo dynamical inconsistency. This can at least be partially removed through 
suitable modifications of the equation which, however, make the method considerably less 
rapid. 

In spite of the above limitations, knowledge of the instability region of the HNC equation 
may allow us to estimate the topology of the region of spinodal decomposition of the fluid. 
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In particular, for the potential defined in Sect. II, with parameters wn/a — 1, wa/cl — 0.2, 
and Ur/Ua = 0.5, it was found that the IL is qualitatively similar to the spinodal line 
calculated through MD calculations^^. More precisely, the density and temperature of the 
low-density critical point estimated through the HNC equation are in satisfactory agreement 
with simulation results, while the density of the second-critical point is overestimated by 
the theory This is not surprising since the theory is an approximate one and becomes 
progressively less accurate as the density increases. However, the ability of the HNC equation 
to give account of the presence of two critical points is, within the well-known limitations 
of the theory, quite remarkable since the potential considered creates a phase diagram that 
is definitely unusual for simple fluids. Thus, studying the modifications of the IL as the 
potential parameters are varied can yield approximate yet useful information on the phase 
behavior of the fluid. In our calculations, functions are evaluated on a grid with M = 2048 
discrete points = mSr, with m — 1, . . . , M and 5r/a — 0.01. 

IV. INSTABILITY LINES FOR THE LARGE REPULSIVE RANGE 

To disentangle the role of each component of the interparticle interaction, the parameters 
of the potential are varied in our investigation one at a time. First we keep the width wr of 
the repulsive shoulder and the width wa of the attractive well fixed and study the behavior 
of the IL, letting the height Ur oi the repulsive shoulder vary. The considered values oHJr 
range from —Ua to oo. When Ur/Ua — —1, the potential consists of a hard core of radius 
a and a square well of width c — a (henceforth called ■potential A) whereas, when Ur ^ oo, 
the potential has a hard core of radius b and a square well of width c — b {potential B). 
When the shoulder height increases, the potential gradually changes from potential A to 
potential B starting from Ur/Ua — —I- In any intermediate configuration, the potential 
has a penetrable finite repulsive shoulder. It is also possible to go from potential A to B 
following a different "path", i.e., by increasing the hard core radius from a to b. In this case, 
however, the potential always consists of an impenetrable hard-core (of varying radius) plus 
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an attractive square well. 

The IL is shown in Fig. 2 at fixed shoulder and well widths {wji/a = 1, wa/o, = 0.2) for 
several values of the shoulder height Ur. In the two limit cases corresponding to potentials 
A and B, the IL exhibits a single maximum corresponding to a phase diagram with a single 
liquid-gas critical point, a well-known behavior for a fiuid of hard spheres with an attractive 
well. The position of the critical points in the p, T plane is considerably different in the two 
cases. The critical point corresponding to potential B is at a lower temperature than that 
corresponding to potential A, due to the weaker attraction, i.e., shorter attractive range c — b 
of potential B with respect to the largest attractive range c — a of potential A. Furthermore, 
the critical density for potential B is smaller than that for potential A and rescales as the 
hard-core volume (a/b)^ of the two potentials. We observe that this rescahng overshadows 
the shift of the critical point toward higher densities due to the decrease of the attraction 
range (e.g., see Appendix A in Ref.^^), unless b/a ~ 1. 

As Ur increases, starting from Ur/Ua — —i (potential A), the IL moves toward lower 
temperatures as a consequence of the overall reduction of the interparticle potential's attrac- 
tive component. At the same time, the IL undergoes a topological change which eventually 
yields a fine with two maxima (Fig. 2b). This peculiar topology of the IL becomes most 
evident for intermediate values of Ur (0.4 < Ur/Ua < 0.6). As Ur increases further, the 
second maximum disappears and again the shape of the IL becomes more and more similar 
to the shape typical of the hard-core square- well potential. Thus, when wr and wa are fixed, 
two maxima are observed in the IL only for a finite range U^""^ < Ur < f/]^*". 

In this range of values, as Ur increases, the density pi of the low-density maximum 
becomes smaller, while that of the other maximum p2 slightly increases. The critical tem- 
peratures Ti and T2 respectively corresponding to these two maxima, decrease — this behavior 
being more evident for the second maximum. These results agree with the behavior found 
with MD simulations for the two critical points (reported in Fig. 9g, 9h of Ref.^^). Thus for 
increasing Ur, the two maxima move away from each other both in density and temperature. 

We now consider a fixed shoulder width {wR/a — 1) and several values of the well width 
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{wa/ol — 0.1, 0.2, 0.3, 0.4, 0.5, 0.6). We calculate the IL for each of them, letting the height 
Ur of the shoulder vary (Figs. 3a, 2, 3b, 4a, 4b and 5a). The values of the shoulder height 
Ur for which two maxima are observed increase with wa- By increasing wa, the second 
maximum is less and less evident and for large values of wa/cl (Figs. 4b and 5a) the second 
maximum is not observed for any Ur. For small values of wa-, the decrease of the attraction 
flattens the curve and the second maximum becomes difficult to observe (Fig. 3a). 

For IL's with two maxima and the same WR/a — 1 and Ur but different wa, both 
maxima move toward higher temperatures for increasing wa-, due to the increased attraction. 
Moreover, by increasing wa, P2 becomes smaller while pi does not vary significantly. This 
behavior agrees with that predicted by MD simulations for the two critical points (shown in 
Fig. 9a,9b of Ref.i^). 

In Fig. 5b we show the behavior of the IL for wr/q = 0.8 and WA/a = 0.2. Comparing 
these results with those shown in Fig. 2, we observe that for fixed wa and Ur, as wr increases, 
Pi and Ti are almost constant and at variance with MD results, but p2 decreases and T2 
increases in agreement with the MD simulations (see Fig. 9d,9e of Ref.^^). 

We next consider a potential with a wider repulsive shoulder {wR/a — 1.5) and several 
values of the well width {wa/ci — 0.5, 1.0, 1.5, 2.5). The behavior of the IL (Figs. 6 and 7) 
for varying Ur is quite similar to that observed in the previous cases. For fixed wr and wa, 
the IL shows only two maxima in a finite range of values of Ur; these values increase with 
Wa and, for large values of wa, the two-maxima topology is not observed regardless of the 
value of Ur (Fig. 7b). The range of values of wa in which we observe two maxima is larger 
with respect to the case wr/q — 1. 

For one particular set of parameters {wr/q — 1.5, wa/o, — 0.5 and Ur/Ua — 0.8), it is 
possible to compare the results obtained using the HNC equation with the phase diagram 
calculated through a theoretical approach based on a thermodynamically consistent integral 
equation^^. Once again it appears evident that the main fiaw of the HNC equation is to 
overestimate the critical density of the second critical point. 

However, a direct comparison of HNC results with those obtained through MD simula- 
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tions can be disappointing. For some of the parameter sets investigated in Fig. 9 of Ref.^^ 
the IL shows only one maximum, while for others the two-maxima topology is barely ob- 
servable. As an example, we show the IL corresponding to the parameters wr/q = 0.5, 
wa/cl — 0.5 with 1.0 < Ur/Ua < 1-7 (Fig. 8). It was not possible to directly analyze the 
value Ur/Ua — 2 (considered in Ref.-*^^) since, in this case, the HNC cannot be solved at high 
densities before any considerable increase of the compressibility can be observed (in general, 
this occurs when the finite repulsion is considerably stronger than the attraction). The re- 
sults obtained at slightly smaller values of Ur show, however, a non-monotonic behavior of 
the IL, suggesting the presence of a liquid-liquid critical point. 



The overall behavior of the IL's is synthesized in Fig. 9 which shows, for different values 
of the shoulder width {wr/g = 0.6, 0.8, 1.0, 1.5), the points in the wa, Ur plane where two 
maxima are found by using the HNC approximation. We observe that both ranges of wa 
and Ur where two maxima are observed increase with wr. 

The general behavior of Ur as a function of wa at constant wr can be rationalized 
by using the modified van der Waals approach (MVDWE) presented in Ref.^^. First we 
approximate the interval of values of C/]^'" <Ur< U'^'"-'-" for each wr and wa in Fig. 9 with 
its middle point U^ = (t/j^""^ + t/]^*")/2 (Fig. 10). Next, we recall from Ref.^^ the relation 
between the potential's parameters and the strength of attraction A, a parameter increasing 
with wa/ci and decreasing with Ur/Ua- In particular, as T — > oo it is 



V. DISCUSSION 



A = Uava - UrVr, 



(1) 



with 



2tt 



[{a + WR + WAf - {a + WR)% 



(2) 



and 



11 



271 

^^ii= Y[(a + w^i?)'-a']- (3) 



The relation Ur/Ua — va/vr — Aj (UaVr) can be rewritten as 



Ur Vsc 



Ua Vsc — Vhc 
where the coefficients 



UaVsc Rsc a Ssc \ a J Vsc \ a J 



(4) 



27r 3 27r 3 Shc a Rhc a 

Vnc^^a ,Vsc= ^{a + WR) ,-— = —— -,— — = — 5) 

3 3 Ssc [a + WRy Rsc a + WR 

are the volumes, and the ratios of the surfaces and radii of the hard core (HC) and the soft 
core (SC), respectively, and all depend only on the parameter wr/g. Hence, at a fixed value 
of wr, the function Ur{wa) in Eq. (4) only has A as an unknown parameter. 

This MVDWE prediction can be verified by using the HNC results. In Fig. 10, Eq. (4) 
is used to fit the values of Uj^{wa) resulting from HNC calculations for different values of 
Wr, with A as the only fitting parameter. As expected from Eq. (4), when wa/o, < 1, the 
leading order in U^{wa) is linear, while when Wa/o, > 1 (corresponding to larger wr), the 
non-linear behavior is evident. Fig. 10 also shows that, by increasing wr, the coefficients of 
the third-degree polynomial in wa decrease as predicted by Eqs. (4) and (5). 

Moreover, the fitting parameter A in Fig. 10 shows a non-monotonic behavior with 
Wr. This is consistent with the MVDWE prediction in Ref.^^ that dA/dwR may have 
different signs, depending on the other parameters. Therefore, Eqs. (4) and (5) give us a fair 
description of how the three parameters Ur, wr, and wa are related to each other when the 
phase diagram has two critical points at positive pressure and finite temperature. However, 
Eqs. (4) and (5) do not help us understand why the phase diagram has an accessible liquid- 
liquid critical point only for limited ranges of wa and Ur, given a value of wr. 

To gain some insight into this point, we observe that if wc plot Eq.(4) with A = and 
no fitting parameters (Fig. 11), we get a rough approximation of the calculated U* that 
becomes fair for the largest wr. This suggests that as a first approximation we can assume 
that A = at least for wr/g > 1, which is consistent with the conclusion of Ref.^^ that in 
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order to have two accessible critical points in the fluid phase, the attractive and repulsive 
part of the potential must compensate, i.e., UaVa — UrVr or ^ ~ 0. 
Hence from Eqs. (2) and (3) we get the approximation 



First we observe that to get an accessible liquid-liquid critical point, Ur/Ua ^ 0{1) is the 
relevant case. Indeed, the case with Ur/Ua ^ 1 at high-enough T and small-enough P 
corresponds to an effective attractive potential with no repulsive shoulder and a hard core 
at a distance of a -|- wr with no liquid-liquid phase transition, or with a liquid-liquid phase 
transition at vanishing T and very high P (see MD results in Fig.s 9h, 9i in Ref.^^). On the 
other hand, for Ur/Ua — 0, Eq. (6) gives wa — 0, leading to a simple hard-core potential 
with no attractive well. Hence we consider the case with Ur/Ua ^ 0{1). 
Next, we observe that for increasing wr {wr — > oo), Eq. (6) goes to 



from which the condition Ur/Ua — ^wa/wr follows for small wa/wr. This relation is 
reasonably satisfied by HNC data for WR/a — 1 {wa/cl ~ 0.4 for Ur/Ua — 0.9) and is better 
approximated for wr/q = 1.5 {wa/cl ~ 0.6 for Ur/Ua — 0.9). We can deduce from these 
considerations that to get a phase diagram with two accessible critical points in the fluid 
region, the three parameters of the potential should be related by the approximate relation 
in Eq. (7) for wr/q 3> 1, which reduces to wa — wr/3 for Ur/Ua — 1- 



The purpose of the present investigation is to understand the role that the different 
components of the interparticle interaction play in the physical mechanism underlying the 
liquid-liquid phase transition in one-component systems. Thus we investigate the phase 
diagram associated with an isotropic pair potential with an attractive well and a repulsive 



Ur ^ (a + WR + wa)^ - (q + wr)^ 
Ua (a + wr)^ - 



(6) 




(7) 



VI. CONCLUSIONS 
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shoulder, described by three parameters, by analyzing for which combinations of the pa- 
rameters do we find a phase diagram with two critical points in the fluid phase. In a flrst 
paper^^ we used MD simulations finding limited ranges of the parameters so that the liquid- 
liquid phase transition was accessible. We also presented a general description based on the 
MVDWE approach, which rationahzed our MD results. 

Since MD simulations are precise but time-demanding, we completed this analysis by 
adopting a different approach with well-known limitations, but one that was extremely fast 
in terms of computational time, consisting of integral equations in the HNC approximation. 
It is important to stress that the drawbacks of the HNC equation are not critical for our 
purposes. Indeed, we found that the theory, though at best only in qualitative agreement 
with MD simulations, correctly reproduces the trend according to which the simulated crit- 
ical points move in the p, T plane as the potential parameters are changed^^. On this basis, 
we use the theoretical results to estimate the phase behavior of our system over a portion 
of the parameter space much wider than that explored by numerical simulations. 

Our findings, both with MD and HNC approach, show that only a limited number of 
combinations of potential parameters can be associated with a phase diagram with two acces- 
sible critical points in the fluid phase. A general conclusion is that the repulsive component 
of the potential must equilibrate the attractive component. By comparing the MVDWE 
predictions with the results presented here, we can quantify the previous general statement 
with the relation 1^41 < UaVr with quantities deflned in Eqs. (l)-(3). 

For wr -C a it is difficult to extract a clear relation between the potential's parameters. 
However, we note that Eq. (4) shows a leading linear relation between Ur/Ua and wa/wr 
for wr <^ a and wa <S a, suggesting that the liquid- liquid phase transition could also be 
found in systems with short repulsive range, if the attractive range is short as well. 

For Wr ^ a the situation is more clear. The MVDWE predictions for A = compare well 
with the HNC results (Fig. 11), leading to the Eq. (6). This equation gives us the intuitive 
understanding that the repulsive and attractive components of the interaction potential 
compensate when the attractive volume, weighted by the attractive energy, is equal to the 
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repulsive volume, weighted by the repulsive energy. Moreover, for large wr, Eq. (6) reduces 
to the simple Eq. (7), whose leading order is Swa/wr. 

In conclusion, we have studied a large number of parameter combinations of a general 
squared isotropic pair potential with an attraction and soft-core repulsion by using the HNC 
equation. We verify the general trend previously observed in MD simulations and extend the 
analysis to a larger number of parameter combinations, with more emphasis on cases with 
wide soft-core repulsion, particularly difficult to simulate by MD. By using the MVDWE 
approach introduced in Ref.^^ to interpret the HNC results, we find that the condition ^4 = 
is well verified for potentials with a large repulsive range and two maxima in the IL. We 
expect that this condition for a liquid-liquid phase coexistence could be generalized to a 
continuous isotropic attractive potential U (r) with a wide soft-core repulsion such as 



where a is the hard-core distance. In particular, this condition brings to Eq. (6) two accessi- 
ble critical points in the fiuid phase between the parameters of squared potentials considered 
here. 




(8) 
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FIGURES 




a b c r 

FIG. 1. General shape of the attractive soft-core potential used in this work, with hard-core 
distance a, soft-core distance b, interaction range c, attractive energy Ua, and repulsive energy Ur. 
We use Ur/Ua, wa = c — b, and wr = b — a as independent parameters. 
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FIG. 2. Instability line of the HNC equation for the potential in Fig. 1 with parameters 
wr/o- = 1, wa/ci = 0.2, and for (from top to bottom in panel a) Ur/Ua = —1, -0.5, 0, 0.2, 
0.4, 0.6, 0.8, 1, oo, and (from top to bottom in panel b) Ur/Ua = 0.3, 0.4, 0.5, 0.6, 0.7, oo. Labels 
A and B marks the curves corresponding to potentials A and B. 
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FIG. 3. Instability line of the HNC equation for the potential in Fig. 1: the parameters of the 
potential are, in panel a, wu/a = 1, wa/ci = 0.1, and (from top to bottom) Ur/Ua = 0.2, 0.3, 0.4, 
0.5; in panel b, wr/g = 1, wa/o, = 0.3, and (from top to bottom) Ur/Ua = 0.4, 0.5, 0.6, 0.7. 
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FIG. 4. Instability line of the HNC equation for the potential in Fig. 1: the parameters of the 
potential are, in panel a, wn/a = 1, wa/o, = 0.4, and (from top to bottom) Ur/Ua = 0.6, 0.7, 0.8, 
0.9, 1; in panel b, wr/q = 1, WA/a = 0.5, and (from top to bottom) Ur/Ua = 0.7, 0.8, 0.85, 0.9, 
1. 
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FIG. 5. Instability line of the HNC equation for the potential in Fig. 1: the parameters of the 
potential are, in panel a, wr/q = 1, wa/o, = 0.6, and (from top to bottom) Ur/Ua = 0.6, 0.8, 0.9, 
1, 1.1, oo; in panel b, wr/g = 0.8, WA/a = 0.2, and (from top to bottom) Ur/Ua = 0.2, 0.3, 0.4, 
0.5. 
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FIG. 6. Instability line of the HNC equation for the potential in Fig. 1: the parameters of the 
potential are, in panel a, wn/a = 1.5, wa/cl = 0.5, and (from top to bottom) Ur/Ua = 0.5, 0.7, 
0.8, 0.9; in panel b, wji/a = 1.5, WA/a = 1, and (from top to bottom) Ur/Ua = 1-5, 1.6, 1.7, 1.8, 
2. 
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FIG. 7. Instability line of the HNC equation for the potential in Fig. 1: the parameters of the 
potential are, in panel a, wr/q = 1.5, wa/o, = 1.5, and (from top to bottom) Ur/Ua = 2, 2.5, 
2.75, 3; in panel b, wr/u = 1.5, WA/a = 2.5, and (from top to bottom) Ur/Ua = 4, 5, 6, 7. 
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FIG. 8. Panel a: Instability line of the HNC equation for the potential in Fig. 1 for wn/a = 0.5, 
wa/cl = 0.5, and (from top to bottom) Ur/Ua = 1, 1-1, 1-2, 1.3, 1.4, 1.5, 1.6, 1.7; panel b: an 
enlarged view at low T. 
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FIG. 9. Symbols mark the combinations of the potential's parameters where the instability line, 
calculated by the HNC approach, shows two maxima, suggesting the presence of two fluid-fluid criti- 
cal points. Sets with wn/a = 0.6 (circles) were investigated for < Ur/Ua < 2 and < wa/cl < 1; 
sets with wr/g = 0.8 (squares) and with wr/q = 1.0 (diamonds) for < Ur/Ua < 1-2 and 
< wa/o- < 0.6; sets with wr/g = 1.5 (triangles) for < Ur/Ua < 6 and < wa/cl < 3. 
Parameters outside these regions have not been investigated. 
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• Wj^/a = 0.6 
□ Wj^/a = 0.8 
w^/a= 1.0 
A Wjj/a=1.5 




FIG. 10. Symbols mark the middle points of intervals of Ur in Fig. 9 for wn/a = 0.6 
(set 1 denoted by circles), wn/a = 0.8 (set 2 denoted by squares), wn/a = 1.0 (set 3 denoted by 
diamonds), wn/a = 1.5 (set 4 denoted by triangles). Error bars represent the interval in Fig. 9. 
Lines are one-parameter fits with Eq. (4) of the sets: for set 2 (squares) the fitting parameter is 
A/ (UaVsc) = 0.31 (dashed line); for set 3 (diamonds) the fitting parameter is A/ (UaVsc) = —0.84 
(dot-dashed line); for set 4 (triangles) the fitting parameter is A/(?7a^5c) = 2.23 (dotted line). 
Since we only have three points for set 1, to avoid a fit with a large indeterminacy on the parameters 
we arbitrarily chose A/ (UaVsc) = 1 to show that the data are consistent with Eq. (4). 
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FIG. 11. Symbols are as in Fig. 10. Lines are Eq.(4) evaluated for ^ = 0: continuous line for 
set 1 (circles), dashed line for set 2 (squares), dot-dashed line for set 3 (diamonds), dotted line for 
set 4 (triangles). 
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